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The Russian farm-fox experiment is an unusually long-running and well-controlled study designed to replicate wolf-to-dog 
domestication. As such, it offers an unprecedented window onto the neural mechanisms governing the evolution of behavior. 
Here we report evolved changes to gray matter morphology resulting from selection for tameness versus aggressive responses 
toward humans in a sample of 30 male fox brains. Contrasting with standing ideas on the effects of domestication on brain size, 
tame foxes did not show reduced brain volume. Rather, gray matter volume in both the tame and aggressive strains was increased 
relative to conventional farm foxes bred without deliberate selection on behavior. Furthermore, tame- and aggressive-enlarged 
regions overlapped substantially, including portions of motor, somatosensory, and prefrontal cortex, amygdala, hippocampus, and 
cerebellum. We also observed differential morphologic covariation across distributed gray matter networks. In one prefrontal-cere- 
bellum network, this covariation differentiated the three populations along the tame-aggressive behavioral axis. Surprisingly, a pre- 
frontal-hypothalamic network differentiated the tame and aggressive foxes together from the conventional strain. These findings 
indicate that selection for opposite behaviors can influence brain morphology in a similar way. 
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Significance Statement 


Domestication represents one of the largest and most rapid evolutionary shifts of life on earth. However, its neural correlates are 
largely unknown. Here we report the neuroanatomical consequences of selective breeding for tameness or aggression in the seminal 
Russian farm-fox experiment. Compared with a population of conventional farm-bred control foxes, tame foxes show neuroana- 


tomical changes in the PFC and hypothalamus, paralleling wolf-to-dog shifts. Surprisingly, though, aggressive foxes also show simi- 
lar changes. Moreover, both strains show increased gray matter volume relative to controls. These results indicate that similar brain 
adaptations can result from selection for opposite behavior, that existing ideas of brain changes in domestication may need revision, 
and that significant neuroanatomical change can evolve very quickly, within the span of <100 generations. 
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Brain volume shifts in the Russian farm-fox experiment. A, Group-average template with anatomic labels and approximate cortical functional localizations from the dog (Evans 


and de Lahunta, 2013). olf, Olfactory; pfc, prefrontal; pmc, premotor; mot, motor; ss, somatosensory; aud, auditory; vis, visual. B, Variation in regional morphology across the dataset. Warmer 
colors represent more variation. Lateral views of 3D surface renderings and coronal cross-sections are shown. C, Strain-wise differences in gray matter, white matter, and total brain volume. 


Error bars indicate SEM. 


Perhaps the most well-known effect of domestication on the 
brain is a reduction in size (Kruska, 2005). Dogs are the oldest 
and perhaps the archetypal domesticate, having split from wolves 
an estimated 10,000-30,000 years ago (Skoglund et al., 2011; 
MacHugh et al., 2017). Past research on the neural correlates of 
wolf-to-dog domestication has implicated PFC and the limbic 
system, particularly the hypothalamo-pituitary-adrenal axis. For 
example, gene expression in the hypothalamus is conserved 
between wolves and coyotes but diverges significantly in dogs 
(Saetre et al., 2004). Similarly, genes showing high differentiation 
between wolves and Chinese native dogs show high expression 
bias for the brain, particularly those expressed in PFC (Li et al., 
2013). In an MRI study of 8 wild carnivore species and 13 
domestic dogs, the allometric scaling of corpus callosum size to 
total brain size was constant across species, except in the rostral 
component, which interconnects prefrontal cortices (Spocter et 
al., 2018). Moreover, the enlargement of one component of car- 
nivore PFC, the prorean gyrus, which has extensive connections 
with other limbic system components (Cavada and Reinoso- 
Suarez, 1985; Markow-Rajkowska and Kosmal, 1987), has been 


implicated in the emergence of complex social behavior in canid 
evolution (Radinsky, 1969). 

The domestication of wolves into dogs is paralleled by 
the well-known and long-running Russian fox experiment 
(Trut et al., 2009). The experiment does not, of course, per- 
fectly recapitulate the “naturally” occurring domestication 
of any species, and indeed domestication trajectories differ 
in important ways across species and can involve selection 
pressures beyond tameness (Zeder, 2012). Nonetheless, the 
Russian fox experiment offers a singularly well-controlled 
window on neurobiological shifts associated with tameness 
(Zeder, 2012). Since 1959, researchers at the Institute of 
Cytology and Genetics in Novosibirsk have been breeding 
conventional farm foxes on the basis of their behavioral 
response to human social contact. The tame strain is 
selected for high social approach behavior toward humans 
and produces dog-like behaviors, such as licking and tail wag- 
ging. In a parallel experiment, the aggressive strain is selected 
for the opposite behavior and reacts with defensive aggression 
when faced with human contact. A third conventional strain 
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Table 1. Trait loadings for PCA of behavior 


Trait 


G7 
A52 
G4 
G1 
G2 
D31 
B12 
A25 
G6 
B25 
B29 
G0 
G3 
B30 
B42 
B13 
B2 
A23 
A40 
G8 
D39 
4 
q 
A31 
GB 
(55 
A32 
G5 
B37 
B48 
B14 
A48 
650 
B9 
A34 
A38 
B8 
A36 
D13 
B39 
Q 
A47 
B47 
B32 
G9 
(204 
D33 
B19 
B3 
D25 
BI 
(25 
6 
C18 
D29 
A37 
C17 
A7 
A29 
B40 
13 
D24 
A2 
D28 


Trait description 


Aggressive sounds 

Did not come to zone 2a 

Follows the hand (aggressive) 

Attack alert 

Pinned ears (aggressive) 

Not on the floor of zone 2 

Not on the floor of zone 2 

Spend in zones 3, 4-6 at least 40 s 

Triangle ears directed back (aggressive) 

Pinned ears (aggressive) 

Spend in zone 3, 4-6 at least 40 s 

Attack 

Trying to bite 

Spend in zones 5 and 6 at least 40 s 

Keeping same posture and place for at least 40 s 
Not on the floor of zone 2a 

Immediately moved back to zone 5 or zones 3, 5, and 6 
Moving back for at least one zone during first 15 s 
Keep same posture and place at least for 40 s 
Animal is present only in zones 3, 5, and 6 
Pinned ears (aggressive) 

Spend more than 30 s in zones 3-6 

First time can touch a fox in zones 5 and 6 

Lie in any zone longer than 30 s 

Animal is in zones 3-6 in the beginning of step C 
Leaning on side or back walls in zones 5 and 6 
Lie in any zone a whole minute 

Narrow ears directed back 

Stays on back or side walls in zones 5 and 6 
Ears are vertical 

Sniffing floor/air 

Sit for at least 20 s 

Tail is up for at least for 3 s 

Sniffing hand from small distance 

Changed place at least once 

Did at least one full circle 

Come to the hand after 10 s 

Changed place at least 2-4 times 

Grooming 

Changed place at least 2-4 times 

Animal is in zones 1 and 2-4 in the beginning of step C 
Tail is up for at least 3 s 

Tail is up for at least 3 s 

Moving on “short legs” 

Moved forward for at least one zone during the step 
Tame sounds (combined) 

Tail is up for at least 3 s 

Loud breathing 

Touch hand for at least 40 s 

Moving forward for at least one zone during first 15 s 
Animal is in zone 2 in the beginning 

Tail wagging 

First time can touch a fox in zones 3 and 4 

Hold hand 

Did at least one full circle 

Changed place at least 5 times 

Rolls on the side, ask to touch belly 

Touch door by foot or scratch 

Came to zones 1 and 2 

Changed place at least 5 times 

Allows to touch back part of the back 

Comes to zones 1 and 2 

Tail wagging 

Changes place at least 5 times 
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PCI 


—0.1476673 
—0.1430533 
—0.1402881 
—0.1399264 
—0.133127 
—0.1304852 
—0.1286657 
—0.1145566 
—0.1084682 
—0.1030547 
—0.0979483 
—0.0811832 
—0.0787731 
—0.0763434 
—0.0746656 
—0.0730084 
—0.0667175 
—0.0635934 
—0.0610294 
—0.0603003 
—0.0589739 
—0.0486048 
—0.0442354 
—0.0413531 
—0.0359193 
—0.0357765 
—0.0322557 
—0.0134525 
—0.0067197 
—0.0051633 
0.00665074 
0.00809613 
0.01117132 
0.01562686 
0.02762113 
0.03153232 
0.03582494 
0.03677646 
0.03813346 
0.0453691 
0.04645646 
0.05655872 
0.05926645 
0.06011035 
0.06502536 
0.06680337 
0.0674776 
0.06814984 
0.06946151 
0.07393203 
0.07482935 
0.07724832 
0.0797726 
0.08141347 
0.084554 
0.08626927 
0.08742998 
0.08991397 
0.09018569 
0.09072512 
0.09318948 
0.0959737 
0.10156213 
0.10330463 


PQ 


0.2057388 
—0.0512895 
0.23863637 
0.2314992 
0.1853784 
—0.0279166 
—0.1034625 
—0.1145345 
0.1639514 
0.12176326 
—0.1150948 
0.20157786 
0.18749222 
—0.1219152 
—0.0280051 
—0.0776614 
—0.1315815 
—0.0798363 
—0.0451188 
—0.2143789 
0.07022877 
—0.1129164 
—0.2005908 
—0.0473826 
—0.0807357 
—0.0621933 
—0.0455808 
—0.1674438 
—0.0593972 
—0.0905863 
—0.0271967 
—0.029441 
0.03423211 
—0.0104944 
0.01619735 
—0.0493083 
—0.0003365 
0.01294125 
—0.0131906 
0.00602985 
0.10165381 
0.00542154 
0.00639307 
0.0280033 
0.13752432 
0.00121025 
—0.0014009 
0.0352152 
0.04423451 
—0.0237983 
0.08558415 
—0.0012107 
—0.0387303 
—0.0092092 
—0.0796958 
0.01591251 
7.A5E-05 
0.0532481 
0.0465907 
—0.0165086 
—0.2048051 
0.01736351 
0.02488945 
—0.0334257 


PC3 


0.05690569 
—0.0320838 
0.07126363 
0.06913457 
0.07517754 
—0.1040545 
0.15847715 
—0.1041976 
0.06852387 
—0.0338697 
0.18041277 
0.02742591 
0.09395528 
—0.040511 
—0.1063443 
0.14535304 
0.16132763 
0.04627051 
—0.0882064 
0.06897475 
—0.0544917 
0.1191827 
0.16847207 
—0.0526588 
0.10769476 
0.09883121 
—0.0514736 
0.07190769 
0.07711811 
0.21988206 
0.12153536 
—0.0223952 
0.11681401 
0.12987889 
0.0878816 
0.12161341 
0.07587808 
0.10609875 
0.00017673 
0.0943546 
—0.0895807 
0.05048287 
0.08498355 
—0.0825452 
0.04538581 
—0.0599009 
0.01885399 
—0.124301 
—0.1249051 
0.10749179 
—0.1237285 
—0.0996494 
0.01207579 
—0.061069 
0.09728685 
0.17384956 
—0.0816944 
—0.0091076 
0.15163891 
0.1744323 
0.00241271 
0.12844973 
—0.0929812 
0.11796783 
(Table continues.) 
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Table 1. Continued 


Trait Trait description 

D6 Touch door by foot or scratch 

B31 Spend in zones 1 and 2 at least 40 s 
B21 Tame ears 

B20 Tail wagging 

D3 Tail wagging 

A24 Spend in zones 1-4 at least 40 s 

A22 Moving forward for at least one zone during first 15 s 
(24 Loud breathing 

A28 Spend in zones 1 and 2 at least 10 s 
A6 Sniffing the door 

(14 Allows to touch back 

AQ Lean on the right wall in zone 2 

D7 Touch door by nose 

D14 Sit in zone 2 and looking at observer 
A10 Sit in zone 2 and looking at observer 
B28 Spend in zones 1-4 at least 40 s 

DI7 Spends in zones 1-4 at least 40 s 
D32 Leaning on right wall in zone 2 

A8 Lean on the door 

A27 Spend in zones 1 and 2 at least 40 s 
(16 Allows to touch head 

AS Touch door with nose 

(8 Lie during a contact for at least 5 s 
B15 Sniffing front door 

C15 Allows to touch nose 

D1 Come to zone 2 during first 5 5 

19 Come to zone 2 in the end of step C 
D2 Spends in zones 1 and 2 at least 30 s 
B7 Trying to poke hand by nose 

D4 Leaning on door 

(29 Come to hand and sniffing in the end of step C 
B11 Come to zones 1 and 2 

C12 Tame ears 

B10 Come to hand and sniffing 
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PC1 PQ PG 
0.10455025 0.01406251 —0.0169301 
0.1063003 0.09768157 —0.1496475 
0.10675677 0.03203084 —0.1522048 
0.10715258 0.03589231 —0.1538206 
0.1081683 —0.0032861 —0.1271162 
0.10833383 0.10252039 0.15880762 
0.11001501 0.0641225 0.15359587 
0.1132731 —0.0010303 —0.1092512 
0.11350046 0.06669391 0.15984224 
0.11575323 0.0628932 0.17880495 
0.11793174 —0.2310268 —0.0218957 
0.11959853 0.07496327 0.07927634 
0.12007958 0.02927169 0.13297053 
0.1224556 0.01267484 —0.0627118 
0.12628208 0.03313335 —0.007816 
0.12730388 0.13902538 —0.1224083 
0.13379193 0.06232127 0.09512188 
0.13470139 0.00703826 0.02006886 
0.1363755 0.06708218 0.07616227 
0.13917388 0.09237863 0.08274502 
0.14020198 —0.2408099 —0.0473998 
0.14037341 0.08793027 0.17366217 
0.14106738 —0.1907798 —0.0817082 
0.14260148 0.10017517 0.05436706 
0.14612813 —0.2334574 —0.0457441 
0.14915423 0.0313106 0.08782056 
0.15108775 0.08114312 —0.0033427 
0.15262612 0.03396873 0.0881678 
0.15360996 0.06452159 —0.0786887 
0.15681275 —0.0184744 0.02205177 
0.15705644 0.02431407 —0.0027936 
0.15826159 0.1000979 0.00695365 
0.1633044 —0.1334769 —0.0702261 
0.16338574 0.06024868 —0.0512387 


Trait loadings for PCA of behavior sorted by PC1. The letter in the trait name indicates a test step (A - observer stands calmly near the closed cage; B - observer opens the cage door, but does not initiate any contact with 
the fox; C - observer attempts to touch the fox; D -observer closes the cage door, then remains calmly near the closed cage) (see Kukekova et al., 2011a for details). 


is kept on the farm but bred without deliberate selection on 
behavior (Statham et al., 2011). Behavioral differences 
between tame and aggressive foxes map to a locus on fox chro- 
mosome 12, which is homologous to the locus implicated in 
wolf-to-dog domestication (Kukekova et al., 2011a, 2012) and 
the whole-genome sequencing of tame, aggressive, and con- 
ventional foxes identifies 103 genomic regions differentiating 
the three populations including 46 regions, which are syntenic 
to canine candidate domestication regions (Kukekova et al., 
2018). Differential gene expression across strains has been 
established in PFC, basal forebrain, hypothalamus, and ante- 
rior pituitary (Kukekova et al., 2011b; Hekman et al., 2018; 
Wang et al., 2018; Rosenfeld et al., 2020). Additionally, tame 
foxes show increased adult neurogenesis in the hippocampus 
in comparison to conventional foxes (Huang et al., 2015). 
However, the brain-wide neuroanatomical consequences of 
selection on behavior in the fox model are as yet unknown. In 
this study, we addressed this question using high-resolution, 
T2-weighted, whole-brain ex vivo neuroimaging in 10 tame, 
10 aggressive, and 10 conventional foxes. These were the same 
individuals used in previous transcriptomic studies (Kukekova et 
al., 2011a, 2012), and their behavior was tested and analyzed as 
previously described (Kukekova et al., 2008). We conducted three 
types of analyses: a comparison of overall gray matter, white mat- 
ter, and total brain volumes; a comparison of differences in re- 
gional gray matter volumes across strains; and a comparison of 


strain-wise differences in anatomic covariation across regions, 
which can identify morphologically coevolving structural net- 
works and link them to individual variation in behavior. 


Materials and Methods 


Brain specimens. In the current study, we examined the brains of 10 
tame, 10 conventional, and 10 aggressive foxes. All foxes were housed 
with littermates and mother from birth until weaning at 6 weeks, then 
housed with littermates for a further 4 weeks, and then singly housed in 
a shed for young foxes with members of the same strain until ~7 months 
of age. After that, all 30 foxes used in this experiment were singly housed 
in the same shed for adult foxes. At all stages of rearing and mainte- 
nance, fox interaction with humans was limited to feeding, veterinary 
care, and testing behavior in home cages at 5.5-6 months of age. All foxes 
were male, sexually naive, and ~1.5 years old at the time of brain extrac- 
tion in August of 2010. Right hemispheres were preserved for gene 
expression studies; we report analyses in left hemispheres here. Left 
hemispheres were formalin-fixed and were maintained in formalin solu- 
tion until scanning in 2016. 

Neuroimaging data acquisition. For imaging, brains were placed in a 
waterproof plastic container, which was packed with polyethylene beads 
for stabilization. The container was then pumped full of Fluorinert FC- 
770 (3M). Fluorinert is a fluorocarbon; it is analogous to a hydrocarbon, 
but with fluorine taking the place of hydrogen. It thus produces no signal 
in (typical) MRI, which is tuned to the resonant frequency of hydrogen 
nuclei. Images were acquired on a 9.4 T/20cm horizontal bore Bruker 
magnet, interfaced to an Avance console, with Paravision 5.1 software 
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(Bruker). A 7.2-cm-diameter volume radio frequency coil was used for 
transmission and reception with a RARE T2 sequence (2 averages, 13 ms 
TE, 2500 ms TR, rare factor 8). Image resolution was 300 um? with a ma- 
trix size of 256 x 100 x 88. 

Image analysis. Image preprocessing was accomplished using the 
FSL software package (Smith et al., 2004; Woolrich et al., 2009; 
Jenkinson et al., 2012). Images underwent bias correction and seg- 
mentation into white matter and gray matter using FAST (Zhang 
et al., 2001). In order to provide a common spatial framework for 
morphometric analysis, an unbiased nonlinear template was built 
from the 10 conventional foxes’ T2-weighted images using the 
ANTS software package (Avants et al., 2009). This template repre- 
sents the group average morphology across the conventional fox 
brain specimens. All subject’s T2-weighted images were nonli- 
nearly aligned to this template. We computed the Jacobian deter- 
minant image of each deformation field; this represents a spatial 
map of where and how much each individual subject’s scan had to 
deform to come into alignment with the template. The Jacobian 
determinant images were then masked with the gray matter seg- 
mentation images to produce representations of each subject’s 
gray matter deviation from the template, and were then smoothed 
with a 0.6 mm (2 voxel) Gaussian kernel. Individual foxes’ gray- 
matter masked, smoothed Jacobian determinant images became 
the input for each of two independent, complementary statistical 
morphometric analyses. The average of all foxes’ masked Jacobian 
images is shown in Figure 1B. 

Voxel-based morphometry (VBM) (Ashburner and Friston, 2000) is 
an inherently hypothesis-based approach that performs GLM at each 
voxel in the image to determine whether morphology is significantly 
related to explanatory variables (e.g., group). VBM analysis was accom- 
plished using FSL’s randomize tool for voxel-wise Monte Carlo permuta- 
tion testing of GLMs, which permutes explanatory variables across cases 
to build up a null distribution, and then tests whether observed associa- 
tions to explanatory variables significantly differ from this random null 
distribution (Winkler et al., 2014). On the other hand, source-based 
morphometry (SBM) (Xu et al., 2009) is a data-driven, model-free 
approach that identifies patterns of significant morphologic correlation 
across subjects. In other words, it determines which regions of the brain 
significantly covary with each other across the entire dataset, while 
remaining agnostic to putative differences across subjects (e.g., group 
differences). Post hoc tests can then be used to determine whether these 
networks show significant associations with variables of interest. This 
was accomplished using the GIFT SBM toolbox for MATLAB (Xu et al., 
2009; http://mialab.mrn.org/software/gift/index.html; for an in-depth 
discussion of the method and its applications, see Gupta et al., 2019). 
Multiple linear regression was then used to assess the relationship between 
factor loadings for covarying gray matter networks and factor loadings for 
a principal components analysis of fox behavioral traits described below. 
Past research has used SBM mainly in the context of identifying neuroana- 
tomical differences between patient and control populations (e.g., schizo- 
phrenia) (Xu et al., 2009; Turner et al., 2012; Rodrigue et al., 2020), but 
several studies have also used SBM in animal studies focused on evolved 
traits: for example, chimpanzee tool use (Hopkins et al., 2019) and eye 
gaze (Hopkins et al., 2020); and breed-specialized skills in domestic dogs 
(Hecht et al., 2019). Covariance in morphology across regions is thought 
to reflect both shared genetic factors and similar patterns of development 
and functional specialization; these structurally covarying networks appear 
to be similar to functionally covarying resting state networks (Alexander- 
Bloch et al., 2013; Evans, 2013; Gupta et al., 2019). Although SBM identi- 
fies networks of structurally covarying regions, it does not assess white 
matter connections linking those regions. 

Assignment of fox behavioral phenotypes. Fox behavior was tested at 
5.5-6months of age in the standard test described by Kukekova et al. 
(2008). Fox behavior during the test was scored from the video records for 
98 recordable observations (Table 1). The matrix, including scores for 
1003 foxes from Kukekova et al. (2011a) and scores for 30 foxes used in 
this study, were subjected to principal component analysis (PCA) in R 
using function: prcomp. The distribution of PC1-PC3 values for 30 foxes 
whose brains are analyzed are shown in Figure 4. 


Hecht et al. e Tameness and Aggression in the Russian Farm-Fox Experiment 


A tame > aggressive B aggressive > tame 


rostral 
composite 
gyrus 


sylvian sulcus 


ppsvy 
ELE 
CPPPS 


D aggressive > conventional 


amygdala hippocampus 


E tame & aggressive both > conventional 


cerebellum hippocampus amygdala premotor cortex & PFC 


Figure 2. _ Differences in regional gray matter volume between strains. A, Tame > aggres- 
sive. B, Aggressive > tame. C, Tame > conventional. D, Aggressive > conventional. E, 
Tame and aggressive both > conventional. Color code: red represents frontal cortex; green 
represents amygdala; magenta represents hippocampus; blue and cyan represents 
cerebellum. 


Results 

The T2-weighted conventional fox brain template and labels 
for some anatomic regions are shown in Figure 1A. A map of 
variation in brain anatomy across all strains is shown in 
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Table 2. Clusters resulting from VBM analysis® 
Cluster Volume 


index (mm?) 
Tame > Aggressive 
1 31.941 
Tame < Aggressive 
2 198.288 
1 
Tame > Conventional 
3 1095.957 
2 303.183 
1 185.22 
Tame < Conventional (none) 
Aggressive > Conventional 
11 2309.283 
10 440.505 
9 104.652 
8 93.717 
7 52.083 
6 22.896 


MAX 
X (vox) 


56 
55 
59 
62 
55 


62 
70 
61 
55 
59 
72 
79 
74 
68 


32 
33 
36 
40 
46 
56 
34 
32 
35 
35 
48 
42 
36 
36 
54 
38 
39 
44 


37 
63 
31 
61 
56 
29 
42 
78 
41 
77 
42 
42 
48 
49 
54 
52 
39 
35 
54 
61 
46 
36 
44 
46 
53 
81 


MAX 
Y (vox) 


162 
160 
140 
132 
150 


189 
186 
196 
193 
175 
180 
139 
141 
154 


222 
223 
222 
226 
207 
210 

98 

85 

73 

86 
100 

96 
129 
151 
122 
141 
145 
133 


222 
190 
223 
198 
202 
217 
130 
123 
121 

96 
140 
143 
149 
141 
149 
145 
148 
129 
125 
122 
141 

90 

86 

77 

79 
132 


(Table continues.) 


MAX 
Z (vox) 


77 
79 
86 
82 
86 


50 
52 
52 
52 
59 


130 
104 
130 


129 
128 


75 
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Table 2. Continued 
Cluster Volume MAX MAX MAX 


index (mm?) X (vox) Y (vox) Z (vox) 
80 139 75 
78 142 65 
5 9.855 52 107 70 
46 103 63 
4 3.402 26 163 58 
30 157 62 
3 3.348 33 74 82 
2 0.567 64 145 59 
1 0.486 42 206 30 


Aggressive < Conventional (none) 


“Each component contains multiple clusters, indicated by Cluster index. Voxel coordinates are for maximum 
value in cluster. 


Figure 1B; a graph of gray matter, white matter, and total 
brain volumes for each strain are shown in Figure 1C. One- 
way ANOVAs revealed a strain-wise difference in total gray 
matter volume (F(227) = 6.855, p=0.004). Post hoc Tukey’s 
HSD tests indicated that both tame and aggressive foxes had 
significantly higher gray matter volume than conventional 
foxes (tame > conventional: t(;g) = 3.370, p = 0.003; aggressive 
> conventional: f(s) = 3.556, p = 0.002), but did not differ sig- 
nificantly from each other (f(13) = 0.010, p = 0.992). Strain-wise 
differences in white matter volume and total brain volume did 
not reach significance (white matter: F(227) = 2.186, p = 0.132; 
total brain volume: Fiy,.97) = 1.504, p = 0.240). 

Voxel-based morphometry identified a number of differences 
between strains (Fig. 2). Carnivore cortical mapping studies have 
been most extensive in cats and ferrets, but relatively less work 
has been done in dogs, which are more closely related to foxes. 
Thus, putative functions of anatomic regions are based on 
known dog anatomy here but should be considered tentative. 
Relative to aggressive foxes, tame foxes show expansion in por- 
tions of the sylvian gyrus and ectosylvian sulcus (temporal 
regions that include auditory cortex and other regions, poten- 
tially higher-order visual or multisensory cortex; Fig. 2A) 
(Kosmal, 2000). In contrast, compared with tame foxes, aggres- 
sive foxes show expansion in portions of the rostral composite 
gyrus and presylvian sulcus (potentially somatosensory and/or 
premotor-prefrontal transition cortex) (Kosmal et al., 1984) and 
in the ectosylvian and sylvian gyri and sylvian sulcus (auditory 
and association cortex in dogs; Fig. 2B) (Kosmal, 2000). 
Surprisingly, relative to the conventional strain, both tame and 
aggressive foxes show expansion in similar regions, including 
portions of the prorean, orbital, frontal, precruciate, and rostral 
composite gyri (prefrontal, premotor, and motor cortex in dogs) 
(Kosmal et al., 1984), amygdala, hippocampus, and cerebellum 
(Fig. 2C,D; overlap shown in Fig. 2E). No voxels showed reduced 
volume in tame or aggressive foxes relative to conventional foxes. 
Anatomical locations, maximum f statistics, number of voxels, 
and volume for each cluster are shown in Table 2. 

Notably, the hypothalamus was one of the regions with the 
highest volumetric variation across the entire dataset (Fig. 1B), 
but direct group-wise comparisons did not identify significant 
volumetric differences between strains (Fig. 2). Past studies on 
the farm-fox experiment have implicated the hypothalamo-pitui- 
tary-adrenal axis generally (for review, see Trut et al., 2009) and 
gene expression in the hypothalamus specifically (Rosenfeld et 
al., 2020). Neuroanatomical consequences of rapid selection on 
behavior can be visible not only in changes to relative volume of 
brain regions, but also in the degree of morphologic correlation 
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Regionally covarying structural brain networks and relationship to strain membership and individual behavior scores. Red-yellow and blue-green components of each network are 


anticorrelated; only significant voxels are shown for each network. Bar graphs represent factor loadings for each strain; positive loadings indicate larger volumes in red-yellow regions (repre- 
sented with red dot), whereas negative loadings indicate larger volumes for blue-green regions (represented with blue dot). Partial regression plots represent significant relationships between 
factor loadings for Networks 2-4 and behavior measurements PC1-PC3 in each fox. A, Network 1 factor loadings did not differentiate between strains and were not significantly related to 
behavior scores. B, Network 2 factor loadings differentiated between tame and aggressive strains and were significantly related to PC1 behavior scores. C, Network 3 factor loadings differenti- 
ated the tame and aggressive strains from the conventional strain and were significantly related to PC2 behavior scores. D, Network 4 factor loadings did not differentiate between strains and 


were marginally related to PC3 behavior scores. 


across brain regions, as we recently documented in domestic dog 
breeds (Hecht et al., 2019). To probe this possibility, we used 
source-based morphometry, a model-free, ICA-based approach 
(Xu et al., 2009), to identify structurally covarying multiregion 
networks across the entire dataset (but white matter connectivity 
is not assessed with this method). This identified four signifi- 
cantly covarying brain networks (Fig. 3; Table 3). 

We then investigated the extent to which these networks were 
related to strain membership and to individual behavior scores. 
Behavior scores were subjected to a PCA analysis, resulting in 
three components, each of which explains 29.0%, 8.0%, and 6.6% 
of variance, respectively. We then used multiple linear regres- 
sion to probe the relationship between factor loadings for 
these three behavior components (Table 1; Fig. 4) and our 
four morphometry components (Table 3). 

Network 1 contained clusters in the thalamus, caudate, NAc, 
cerebellum, and other regions (Fig. 3A). Factor loadings were in 


opposite directions for the tame and aggressive strains but were 
centered near zero with wide variance in the conventional strain; 
however, strain-wise differences did not reach significance 
(Fi2,27) = 0.475, p=0.627). Multiple linear regression indicated 
that factor loadings for Network 1 accounted for 11.6% of the 
variance in behavior PCA scores, but the overall regression 
model did not reach significance (F(3,26) = 1.137, p=0.352), nor 
did any of the partial correlations between brain and behavior 
factor loadings. 

Network 2 was mainly comprised of one cluster that covered 
most of the hypothalamus, plus several other clusters scattered 
throughout the cerebellum. Factor loadings strongly differenti- 
ated the tame and aggressive strains from each other and were 
centered near zero for the conventional strain. Strain-wise differ- 
ences were significant (F227) = 4.495, p = 0.021; Fig. 3B). Post hoc 
tests using Tukey correction for multiple comparisons confirmed 
that factor loadings significantly differentiated the tame and 
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Table 3. Clusters resulting from SBM analysis 


Network 1 (negative/blue) 


Network 1 (positive/orange) 


Cluster 
index 


18 


17 


15 


14 


me Nw en 


= 


4 


PNW BNW ~ C 


Volume 


(mm’) 


137.997 


125.712 


48.573 
39.312 


17.523 


8.451 
6.129 
5.265 
5.238 
5.184 
2.349 
1.755 
1.188 


1.053 
0.837 
0.567 
0.459 
0.432 


383.94 


94.338 


41.499 
24.516 


17.604 
16.659 


6.831 
4.617 
3.996 
3.402 
2.052 
0.864 
0.162 
0.027 


MAX 
X (vox) 


MAX 
Y (vox) 


114 
116 

95 
126 
150 
135 
145 
149 
160 
158 
152 

80 
218 
210 
204 
219 
208 
202 
213 
215 

87 
118 

86 
122 
156 
195 
192 
190 
209 
155 
233 
107 

99 


51 

55 

53 

56 

49 
102 
234 
229 
225 
238 
234 
214 

73 

61 
105 
108 

99 
118 
113 
107 
201 
178 
101 
186 
134 
142 
157 
151 
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Table 3. Continued 


MAX Cluster Volume 
Z (vox) index (mm?) 
Network 2 (negative/blue) 
120 9 153.954 
114 
122 
122 8 40.635 
123 7 39.501 
68 
63 
70 
61 
56 
72 6 33.966 
51 
53 5 25.596 
34 
44 
87 4 12.15 
97 
97 
70 3 9.18 
19 2 0.378 
48 1 0.054 
86 Network 2 (positive/orange) 
68 11 528.201 
50 
53 
57 
61 
62 
85 10 48.168 
51 
30 
56 9 8.289 
69 8 7.128 
7 5.373 
89 6 1.89 
83 5 1.674 
70 4 0.945 
75 3 0.648 
68 2 0.027 
78 1 0.027 
57 Network 3 (negative/blue) 
49 4 4.428 
73 3 2.241 
59 
69 2 0.729 
44 1 0.324 
42 Network 3 (positive/orange) 
55 4 898.344 
60 
54 
44 
103 
100 
100 3 11.205 
94 2 3.159 
36 1 0.513 
77 Network 4 (negative/blue) 
102 14 14.31 
51 
57 13 12.285 
45 
40 12 9.531 


(Table continues.) 


53 


MAX 
Y (vox) 


55 
88 
97 

158 

1733 

156 

158 

168 

133 

131 

161 

154 
97 
85 
82 

232 

237 

236 
44 
99 

112 


76 
84 
98 
55 
95 
56 
90 
92 
90 
145 
75 
118 
135 
223 
114 
216 
218 
63 


90 
47 
43 
95 
89 


225 
217 
212 
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203 
139 
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103 
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163 
79 


MAX 
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93 
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Table 3. Continued 


Cluster Volume MAX MAX MAX 
index (mm?) X (vox) Y (vox) Z (vox) 

52 84 49 

11 9.018 76 121 55 

10 4.239 49 98 64 

9 2.106 32 238 62 

8 1.377 45 64 66 

7 1.188 32 233 75 

6 0.918 39 107 45 

5 0.783 60 87 67 

4 0.621 50 100 82 

3 0.216 28 125 74 

2 0.081 54 74 52 

1 0.027 74 107 53 

Network 4 (positive/orange) 

a3 452.682 24 156 38 

21 139 36 

26 130 41 

26 143 59 

21 127 35 

21 152 52 

11 184.788 17 97 53 

Fy 114 71 

20 90 74 

16 126 60 

5 112 65 

16 111 55 

10 98.685 26 46 68 

25 46 51 

24 46 49 

9 26.109 aw 117 130 

26 112 129 

28 107 128 

zi 133 131 

ai 141 130 

8 10.773 24 65 78 

7 10.746 23 67 50 

6 1,539 36 94 122 

5 0.486 32 135 92 

4 0.27 21 63 38 

3 0.162 56 215 63 

2 0.081 21 66 94 

1 0.027 35 134 85 


aggressive strains from each other (p=0.016). Loading coeffi- 
cients for Network 2 accounted for 31.0% of the variance in 
behavior scores. The overall regression model reached signifi- 
cance (F(3.26) = 3.908, p=0.020). In examining partial correla- 
tions with individual behavior categories, we applied the 
Benjamini-Hochberg method for FDR correction. This revealed 
a significant partial correlation with PC1 behavior scores (8 = 
0.166, t=2.964, p =0.006, which exceeded the critical value of 
0.017). Generally, behavioral traits with the highest positive load- 
ings for PC1 describe tame behavior and proximity to a human 
approacher, whereas traits with the most negative loadings 
describe aggressive behavior and greater distance from a human 
approacher (Table 1; Fig. 4). Examination of the partial correla- 
tion scatter plot reveals that tame foxes cluster together with 
higher scores for both behavior and neuroanatomy, whereas 
aggressive and conventional foxes cluster together on the lower 
end of both axes (Fig. 33). 

Network 3 contained two smaller, discrete clusters in the 
hypothalamus, plus a large cluster covering much of prefrontal 
and premotor cortex, including portions of the prorean, orbital, 
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frontal, precruciate, and rostral composite gyri. For this network, 
factor loadings were similar for the tame and aggressive strains, 
and both were strongly opposite to the conventional foxes; the 
strain-wise difference was significant (F(2,27) = 14.795, p < 0.001; 
Fig. 3C). Post hoc Tukey tests revealed that this network signifi- 
cantly differentiated both the tame strain from the conventional 
strain (p = 0.001) and the aggressive strain from the conventional 
strain (p=0.000). Regression analyses revealed that behavior 
scores accounted for 35.3% of the variance in morphometry 
loading coefficients. The overall regression model reached signif- 
icance (F(3,26) = 4.720, p = 0.009). Examination of partial correla- 
tions using Benjamini-Hochberg FDR correction revealed that 
there was a significant partial correlation with PC2 (8 = 0.483, 
t=3.251, p=0.003, which exceeded the critical value of 0.017). 
For PC2, traits with the highest positive loadings describe active 
aggressive response and fox position in the front part of the cage, 
whereas traits with highest negative loadings describe a passive 
tame response: tolerance of human tactile contact and fox posi- 
tion in the back of the cage (Table 1; Fig. 4). Examination of the 
partial correlation scatter plot reveals that aggressive foxes 
tended to score high for both the behavior and neuroanatomical 
measures, whereas conventional foxes tended to score low on 
both, and tame foxes were tightly clustered in the intermediate 
range (Fig. 3C). 

Network 4 contained a large cluster that spanned regions of 
the thalamus, most of the hypothalamus, and portions of the 
NAc/ventral forebrain and caudate (Fig. 3D). There were also 
several additional clusters located in cortex and cerebellum. 
Factor loadings were again opposite for the tame and aggressive 
strains, and centered near zero for the conventional strain; 
strain-wise differences did not reach significance (F227) = 1.198, 
p=0.312). Loading coefficients for Network 4 accounted for 
14.6% of the variance in behavior scores. The overall regression 
model did not reach significance (Fi3.6) = 1.413, p=0.242). 
None of the partial correlations between morphometry and 
behavior scores reached significance, but there was a marginally 
positive relationship with PC3 (8 = 0.492, t=1.780, p= 0.087, 
which did not exceed the critical value of 0.017). In PC3, traits 
with the highest positive loadings describe neutral exploratory 
behavior (e.g., “ears are vertical”), whereas traits with the most 
negative loadings describe prosocial greeting behavior and fox 
position in the front part of the cage (see Table 1). Tame, aggres- 
sive, and conventional foxes showed marked overlap in both 
behavior and neural scores in this analysis (Fig. 3D). 


Discussion 

The Russian farm-fox experiment is perhaps the longest-run- 
ning, best-controlled, and most well-known artificial selection 
study bearing on the evolution of mammalian behavior. As such, 
it enables a uniquely powerful window on the neural mecha- 
nisms governing behavioral adaptation. This study used high-re- 
solution MRI to examine the brains of these foxes. We found 
that selection on social behavior has altered the anatomy of dis- 
tributed gray matter networks, which included, among other 
regions, PFC, hippocampus, amygdala, caudate, NAc, cerebellum, 
and hypothalamus. These regions have also been implicated in 
past studies in these foxes (Kukekova et al., 2011b; Huang et al., 
2015; Hekman et al., 2018; Wang et al., 2018; Rosenfeld et al., 
2020) and in past studies on wolf-to-dog domestication (Nikulina, 
1991; Saetre et al., 2004; Natt et al., 2012; Li et al., 2013; Ruan and 
Zhang, 2016; Spocter et al., 2018; Oshchepkov et al., 2019). 
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Figure 4. Population distributions for the first three PCs (PC1-PC3) of fox behavior. Aggr, “Aggressive” population; Conv, “conventional farm-bred population’; Tame, “tame” population. 


Black dots represent individual data points. Horizontal bars within each box represent the population median. The whiskers represent the range of data up to 1.5 times the interquartile range. 


Interestingly, some of these regions are also affected by domestica- 
tion in rabbits (Brusini et al., 2018). 

Notably, portions of fox PFC appeared in several of our 
results (Figs. 2E, 3A, 3C). How can the functional relevance of 
this region be interpreted? Although some of these prefrontal 
results are located on the dorsolateral surface of the brain (i.e., 
the prorean gyrus), this region is likely not homologous to the 
granular dorsolateral PFC of humans and macaque monkeys, 
as that region is thought to be unique to primates (Preuss, 
1995; Passingham and Wise, 2012). Rather, the connectivity 
and cytoarchitecture of carnivore PFC (Narkiewicz and 
Brutkowski, 1967; Kosmal and Dabrowska, 1980; Kosmal, 
1981a,b; Stepniewska and Kosmal, 1986; Markow-Rajkowska 
and Kosmal, 1987; Rajkowska and Kosmal, 1988) are more 
similar to that of the dysgranular and agranular portions of 
primate orbitofrontal and ventromedial PFC, regions that 
function to integrate external multisensory and internal vis- 
ceromotor information with limbic reward and threat signals 
to compare the value of potential behavioral choices 
(Carmichael and Price, 1996; Rudebeck and Murray, 2011). 
Putatively, this circuitry has been altered in tame foxes to 
bias behavioral decisions toward the reward value of social 
contact with humans, and in aggressive foxes toward the op- 
posite. This interpretation fits generally with conceptualiza- 
tions of partially shared circuitry between social approach 
and social avoidance processing in humans and other species 
(Aupperle and Paulus, 2010). 

However, the current study also produced some findings 
which were unexpected and suggest revision of existing thinking 
about domestication. One of these was that total gray matter vol- 
ume in both the tame strain and the aggressive strain is increased 
compared with the conventional strain. These findings are in 
contrast to a number of prior studies which have reported that 
domestication reduces brain size in diverse species including 
Atlantic cod (Mayer et al., 2011), guppies (Burns et al., 2009), 
rainbow trout (Marchetti and Nevitt, 2003), mallard ducks 
(Ebinger, 1995), rats, mice, gerbils, guinea pigs, rabbits, pigs, 
sheep, llamas, horses, ferrets, cats, and dogs (for review, see 
Kruska, 2005). Brain size should be interpreted in the context of 
body size, given that the two covary. While body size measure- 
ments were not available for the foxes in our study, Huang et al. 
(2015) report a 15.9% reduction in body weight in tame com- 
pared with conventional foxes, suggesting that, if anything, dif- 
ferences in brain:body size ratios are more pronounced than the 
brain measurements alone would indicate. What might cause 
this increase in brain size? One possible explanation might 


involve neuroplastic changes resulting from differential lifetime 
experience. The link between naturalistic, enriched environments 
and larger, more complex brains has been noted for decades (van 
Praag et al., 2000; Lambert et al., 2019). Tame, aggressive, and con- 
ventional foxes are all housed in identical conditions and undergo 
identical treatment with minimally necessary human interaction. 
Nonetheless, tame foxes, because of their innate predisposition to- 
ward prosocial interaction with humans, may effectively experi- 
ence this environment to be more naturalistic and enriched. 
However, this potential interpretation does not explain 
why the aggressive strain also shows increased gray matter vol- 
ume relative to controls. Perhaps their constant drive to avoid 
human contact functions as a sort of enrichment, but we instead 
favor an alternative explanation. Both the tame and aggressive 
strains have been subject to intense, sustained selection on 
behavior, while the conventional strain undergoes no such inten- 
tional selection. Thus, it is possible that fast evolution of behav- 
ior, at least initially, may generally proceed via increases in gray 
matter. Several potential mechanisms might underlie this effect. 
First, because of developmental linkages, selection pressure for 
increased size in one region may “drag along” enlargement in 
others to some extent (Finlay et al., 2001), although it seems 
apparent that the current dataset also reflects at least some degree 
of modular evolution (i.e., focused selection pressure on specific 
regions and networks) (Striedter, 2020). Second, selection pres- 
sure for expansion of a particular brain region might favor 
genomic variants leading to changes in expression of genes 
with potential pleiotropic effects (e.g., transcription factors). 
Consequently, these gene expression changes could produce 
expansion not only in that one region but also in others. 
Normally, trade-offs with metabolic and life history constraints 
(Isler and van Schaik, 2009) might favor further fine-tuning of the 
activity of these genes so that the effect occurs more specifically 
only in the “targeted” brain region, but because captivity reduces 
the pressure for optimization, this fine-tuning might fail to occur. 
Third, because of the complex interdependencies between brain 
systems, adaptive solutions that modify one existing system may 
stand a high likelihood of producing deleterious effects in others 
unless additional compensatory adaptations also occur (De Vries, 
2004). In the context of extreme selection pressure, these compen- 
satory changes might be more likely to occur via the addition of 
new neural material rather than via volumetrically net-zero altera- 
tions to existing circuits, which putatively would require more 
fine-tuned genetic changes. In the wild, there would be constant 
pressure against these “easy, wasteful” solutions, but in the context 
of the farm-fox experiment, these constraints may have been 
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lessened or removed, as these foxes are not bred for increased 
meat yields, more rapid maturation, or reduced feeding costs. In 
general, it is important to note that strong selection for a single 
dimension of behavior in captivity, as occurred in the farm-fox 
experiment, is expected to differ in important ways from selection 
in the wild, which occurs on many behaviors at once. Because of 
this, experimentally applied selection pressure might be more 
likely to produce “nonoptimal” brain network changes, in contrast 
to in the wild, where such shifts might produce detrimental effects 
on other behaviors, and, importantly, on animal survival. There 
are multiple potential cellular-level causes for increased gray mat- 
ter, including increased neuron count, increased neuron size, and/ 
or dendritic changes, such as increased arborization or spine den- 
sity (Zatorre et al., 2012; Keifer et al, 2015). Future research, 
including single-cell transcriptomics and histologic work, will be 
required to differentiate among these possibilities. 

Notably, a recent skeletal morphology study did not find a 
significant difference in endocranial volume between tame and 
conventional foxes (Kistner et al., 2021). The foxes that became the 
progenitor population for Belyaev’s study had existed in a farmed 
state in Russia for ~50 years (Kochergin, 1936; Vahrameyev and 
Belyaev, 1948); this original Russian farmed population was itself 
drawn largely from Eastern Canadian foxes (Statham et al., 2011; 
Rando et al., 2017) bred at fur farms on Prince Edward Island 
(Forester and Forester, 1982). Kistner et al. (2021) compared endo- 
cranial volumes of modern Russian farm-fox skulls to skulls from 
wild foxes, which had been collected in Canada east of Quebec 
between 1884 and 1952, with 70% collected between 1894 and 1900. 
This revealed that endocranial volume in the modern, conventional 
Russian farm-raised foxes (ie. those bred without selection on 
behavior) was significantly reduced compared with the archived 
Canadian wild skulls. However, no significant difference was 
observed between modern tame and conventional strains, but nota- 
bly, like the current study, Kistner et al. (2021) were unable to incor- 
porate body size measurements into their analyses. We also do not 
know whether the skull morphology of conventional foxes started 
to show differences at the time when selection of conventional foxes 
for tameness in Russia began (1959) or whether it became pro- 
nounced later when the conventional population was evolving in 
parallel with the tame population. Together with these findings, the 
current study hints that the brain volume reduction associated with 
domestication might occur relatively quickly, before the onset of 
intentional selection on behavior. Further research is required to es- 
tablish whether this is true. 

A second surprising result from this study was that both tame 
and aggressive foxes showed enlargement in substantially overlap- 
ping gray matter regions, including PFC, amygdala, hippocampus, 
and cerebellum (Fig. 2E). In addition to these volumetric effects, 
we also observed strain-wise differences in the degree of morpho- 
logic covariation across distributed, multiregion networks. These 
latter measurements revealed links with individual variation in 
behavior, including in brain regions that did not show volumetric 
differences, notably, the hypothalamus. Additionally, Network 3 
consisted primarily of the hypothalamus and PFC, two regions 
strongly implicated in both fox and dog domestication (Nikulina, 
1991; Saetre et al., 2004; Kukekova et al., 2011b; Natt et al., 2012; 
Li et al., 2013; Ruan and Zhang, 2016; Spocter et al., 2018; Wang 
et al., 2018; Oshchepkov et al., 2019; Rosenfeld et al., 2020). In this 
network, factor loadings did not differentiate the tame from 
aggressive strains; rather, the selectively bred strains together were 
differentiated from the conventional strain (Fig. 3C). Together, 
these results indicate that selection for opposite behavioral 
responses (docility vs aggression) can produce similar evolved 
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changes in the brain. This has important implications for attempts 
to evaluate the hypothesis that humans are ourselves self-domesti- 
cated (Sanchez-Villagra and van Schaik, 2019), given that similar 
neuroanatomical patterns of change could now be interpreted 
to support either selection for increased or decreased aggression in 
our lineage. Notably, Trut et al. (2017) reported some elements of 
the “domestication syndrome” in aggressive foxes, including 
altered timing of annual reproductive activity, facial foreshorten- 
ing, and facial “star” marking; they thus propose that these traits 
might result pleiotropically from selection on either tameness or 
aggression (i.e., as a result of selection on social behavior gener- 
ally). This calls for additional research on the neural correlates of 
evolved differences in both docility and aggression at the cellular 
and genomic levels. 

The Russian farm-fox experiment represents a uniquely well- 
controlled opportunity to study the effects of specific, sustained 
selection on behavior. Thus, apart from questions of domestica- 
tion, an additional implication of these results concerns brain 
evolution on a more general level. We found that intense selec- 
tion on behavior can produce gross changes in distributed brain 
morphology extremely rapidly (within the span of well under 
100 generations). This suggests that the brains of other animals 
on this planet, including Homo sapiens, may have undergone 
similarly precipitous morphologic shifts any time steep selection 
on behavior was experienced. 

“Within each network, the negative components (shown as 
blue in Fig. 3) and positive components (shown as red in Fig. 3) 
are morphologically anticorrelated. Each component contains 
multiple clusters, indicated by Cluster index. Voxel coordinates 
are for maximum value in cluster. 
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